***Figure 6.1, Caloric Potential Cell***

version 15.1

clear
use "worms_replication.dta"
set more off

***generate time counter***
gen time = year - 324

***set caloric potential as wealth measure
gen wealth =  calories
label var wealth "Wealth"
	
*Create half Century Counters*
gen half_century = 50*floor(year/50)
label var half_century "Half Century"
quietly forval i = 3/18{
	gen HC`i'0 = 0
	gen HC`i'5 = 0
	replace HC`i'0 = 1 if year >= `i'*100 -50 & year <`i'*100
	replace HC`i'5 = 1 if year >= `i'*100 & year <`i'*100 + 50
	}

*Create kingdom by half-century identifiers***
gen kingdom_hc =  half_century + .0001* kingdom_id 
gen grid400_hc = half_century + .001*grid_group_400

	
/* Collapse data to diocese-bishop observations */
sort diocese year
drop if byearin != year


	postutil clear
	postfile Results coefficient upper lower startyear using "Results.dta", replace

	***Generate the placebo treatments and run the regression for each year, saving coefficient and standard error***
	preserve
	 forval  i = 800/1275{
		gen treat_period = 0
		//generate the placebo treatment periods for regions covered by Worms
		replace treat_period = 1 if year >= `i' + 15 & year <=`i' + 202
		replace treat_period = 0 if year >=`i' & year < `i' + 15 & Burgundy == 1
		replace treat_period = 0 if year >=`i' & year < `i' + 15 & country == "France" & kingdom_id == 14
		//generate the placebo treatment periods for regions covered by London and Paris
		replace treat_period = 1 if year >=`i' & year <= `i' + 202 & country == "England_and_Wales"
		replace treat_period = 1 if year >=`i' & year <= `i' + 202 & country == "France"
		replace treat_period = 1 if year >=`i' & year <= `i' + 202 & diocese == "Tournai"
		replace treat_period = 0 if country == "Ireland"
		replace treat_period  = 1 if country == "Ireland" & year >=`i' + 64
		gen treat = treat_period * covered
		gen interaction = wealth * treat
		quietly reghdfe  religious_bishop   interaction treat crusade1 crusade2 crusade3 crusade4, a(iddiocese kingdom_hc) vce(cluster iddiocese)
		post Results (_b[interaction])  (_b[interaction] + _se[interaction]*invttail(e(df_r),.025))  (_b[interaction] - _se[interaction]*invttail(e(df_r),.025)) (`i')
		drop treat_period treat interaction
			}
	
	***Draw the figure using the stored coefficients and standard errors***
	postclose Results
	clear
	use "Results.dta"
	label var startyear "Start Year for Treatment Period"
	label var coefficient "Estimated Coefficient for Interaction Term"
	label var upper "Upper Confidence Interval"
	label var lower "Lower Confidence Interval"
	graph twoway line coefficient upper lower startyear, xlabel(800(50)1275)  yline(0, lpattern(dot) lcolor(black))  xline(1107, lpattern(dot) lcolor(black)) mcolor(black black black) lcolor(black gs8 gs8) lwidth(thick medium medium)  graphregion(color(white))  ytitle(Coefficient on Interaction) ylabel(,nogrid) bgcolor(white) legend(off) title("Caloric Potential")
	postutil clear	
	restore 

	
